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Intelligent  energy  management  is  a  cost-effective  key  path  to  realize  efficient  automotive  drive  trains  [R. 
O’Hayre,  S.W.  Cha,  W.  Colella,  F.B.  Prinz.  Fuel  Cell  Fundamentals,  John  Wiley  &  Sons,  Hoboken,  2006],  To 
develop  operating  strategy  in  fuel  cell  drive  trains,  precise  and  computational  efficient  models  of  all  system 
components,  especially  the  fuel  cell  stack,  are  needed.  Should  these  models  further  be  used  in  diagnostic 
or  control  applications,  then  some  major  requirements  must  be  fulfilled.  First,  the  model  must  predict 
the  mean  fuel  cell  voltage  very  precisely  in  all  possible  operating  conditions,  even  during  transients.  The 
model  output  should  be  as  smooth  as  possible  to  support  best  efficient  optimization  strategies  of  the 
complete  system.  At  least,  the  model  must  be  computational  efficient.  For  most  applications,  a  difference 
between  real  fuel  cell  voltage  and  model  output  of  less  than  10  mV  and  1000  calculations  per  second 
will  be  sufficient.  In  general,  empirical  models  based  on  system  identification  offer  a  better  accuracy 
and  consume  less  calculation  resources  than  detailed  models  derived  from  theoretical  considerations  [J. 
Larminie,  A.  Dicks.  Fuel  Cell  Systems  Explained,  John  Wiley  &  Sons,  West  Sussex,  2003  ].  In  this  contribution, 
the  dynamic  behaviour  of  the  mean  cell  voltage  of  a  polymer-electrolyte-membrane  fuel  cell  (PEMFC) 
stack  due  to  variations  in  humidity  of  cell’s  reactant  gases  is  investigated.  The  validity  of  the  overall  model 
structure,  a  so-called  general  Hammerstein  model  (or  Uryson  model),  was  introduced  recently  in  [M. 
Meiler,  0.  Schmid,  M.  Schudy,  E.P.  Hofer.  Dynamic  fuel  cell  stack  model  for  real-time  simulation  based  on 
system  identification,  J.  Power  Sources  176  (2007)523-528],  Fuel  cell  mean  voltage  is  calculated  as  the  sum 
of  a  stationary  and  a  dynamic  voltage  component.  The  stationary  component  of  cell  voltage  is  represented 
by  a  lookup-table  and  the  dynamic  voltage  by  a  parallel  placed,  nonlinear  transfer  function.  A  suitable 
experimental  setup  to  apply  fast  variations  of  gas  humidity  is  introduced  and  is  used  to  investigate  a  10 
cell  PEMFC  stack  under  various  operation  conditions.  Using  methods  like  stepwise  multiple-regression  a 
good  mathematical  description  with  reduced  free  parameters  is  achieved. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cell  vehicles  are  a  promising  alternative  to  vehicles  pow¬ 
ered  by  internal  combustion  engines  because  of  their  low  noise  and 
zero-emission.  The  core  piece  of  a  fuel  cell  power  train  is  a  stack  of 
polymer-electrolyte-membrane  fuel  cells  (PEMFCs)  supplied  with 
hydrogen  and  oxygen. 

One  of  the  major  issues  in  the  operation  of  PEMFC  stacks  is 
a  satisfying  humidification,  since  the  humidification  level  of  the 
stack  plays  a  decisive  role  for  performance.  First  of  all,  the  proton 
conductivity  of  the  membrane  is  proportional  to  the  water  con¬ 
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tent  [1],  The  membrane  is  a  derivation  of  polytetrafluoroethylene 
(PTFE)  with  an  additional  compound  of  H+-S03-.  In  the  presence 
of  water,  this  bond  is  just  weak  and  channels  for  motile  protons 
will  be  created  [2],  The  protons  drag  water  molecules  from  anode 
to  cathode  by  an  electro-osmotic  effect  and  in  consequence,  the 
anode  contains  less  water  than  the  cathode,  where  furthermore 
water  is  produced.  However,  also  the  cathode  will  dry  out  at  high 
temperatures  or  high  flow  rates.  To  avoid  this  drying,  to  balance 
the  different  humidification  levels  of  the  electrodes,  and  to  mini¬ 
mize  the  overvoltage  caused  by  membrane  resistance,  the  gas  must 
be  humidified  before  entering  the  stack.  But  even  an  excess  of 
humidification  can  decrease  the  stack  performance  if  water  floods 
the  electrodes,  blocks  the  pores  in  the  gas  diffusion  layer,  and  so 
decrease  the  mass  transport  significantly.  Therefore,  a  controlled 
water  balance  is  an  important  task  for  the  operation  of  a  PEMFC 


f  al. /Journal  of  Power  Sources  190  (2009 )  56-63 


In  transportation  applications,  the  stack  voltage  is  dependent 
of  several  surrounding  conditions  and  80%  of  operating  situations 
are  dynamic  [3],  The  development  of  new  operating  strategies  is 
therefore  a  challenging  task.  In  particular  the  influence  of  fuel  cell 
gas  humidity,  a  high  relation  to  dynamic  behaviour  of  stack  voltage 
can  be  found. 

So  far,  only  models  for  the  static  behaviour  of  stack  voltage  in 
relation  to  gas  flow  humidity  are  known  in  literature  [4,5].  In  this 
work,  a  valid  model  for  the  simulation  of  the  dynamic  behaviour 
of  stack  voltage  by  variation  of  fuel  cell  gas  humidity  is  introduced 
and  developed  for  the  first  time.  The  aim  of  this  work  is  to  provide 
a  simulation  tool  which: 

•  simulates  the  stack  voltage  with  a  precise  prediction, 

•  is  valid  in  power  train  relevant  operating  area, 

•  uses  an  easy  and  reliable  identification, 

•  is  suitable  for  simulations  in  real-time. 

2.  Theoretical  background 

2.J.  Gas  humidification 

The  saturation  pressure  E  of  steam  can  be  approximated  with 
the  actual  steam  inlet  temperature  T  (in  °C)  by  using  the  Magnus 
equation  for  water  [6]  as  given  in  the  following  equation: 

E  =  610.78  Pa  x  10f7-5T«237-3+r»  (1) 

Using  this  saturation  pressure,  the  steam  pressure  e  as  func¬ 
tion  of  the  desired  relative  humidity  rh  can  be  calculated  using  the 
following  equation: 

e  =  rh£  (2) 

The  mixing  ratio  mr  is  defined  by  the  ratio  of  mass  flow  of  steam 
mHjO  t0  mass  fl°w  °f  dry  gas  mgasdry  or  with  the  molar  masses 
Mgas.dry.  Mh2o.  and  the  absolute  pressure  of  the  gas  p,  respectively: 

mr  -  _  mh2o  e  ^ 

nigas.dry  Mgas.dry  p  —  e 

Thus,  the  required  steam  mass  flow  can  be  calculated  by  the 
following  equation: 

mH2o  =  mrmgas>dry  (4) 


probability  of  the  null  hypothesis  0  is  smaller  than  entrance  level 
0in  (0  <  0in)  the  null  hypothesis  has  to  be  rejected  and  the  tested 
term  is  included  into  the  model.  As  result,  a  reduced  linear  model 
with  only  the  most  significant  parameters  is  available  to  predict  the 
output  variable.  If  there  are  more  model  terms  <j>  than  experimental 
data  0  available,  not  included  terms  are  added  to  model  individu¬ 
ally  and  their  0-values  are  calculated  successively.  The  term  with 
the  smallest  0-values  is  added  to  the  model.  This  is  repeated  until 
no  excluded  terms  have  a  0-values  below  0out- 

Following,  the  used  equations  during  stepwise  regression  are 
described  briefly.  First,  the  optimal  parameters  |  of  the  linear 
model: 

y  =  X|  (5) 

have  to  be  estimated.  X  is  a  0-by-0  matrix  of  the  input  variables  and 
Y  the  vector  of  the  estimated  output  variable  K d.  The  estimation  of 
the  model  quality  is  possible  by  the  square  error  amount  J  between 
the  estimated  and  measured  output  variable: 

d>  <p 

J  =  ^2(Yd-Yd)2  =  ^(  Yd  -  Xdt)2  (6) 


2.2.  Stepwise  regression  analysis 

A  structured  process  to  determine  an  empirical  model  is  the 
stepwise  regression  method  [7],  The  aim  of  this  method  is  to  gain 
a  linear  model  equation  for  the  target  variable,  which  is  as  sim¬ 
ple  as  possible  but  still  estimates  the  output  variable  with  a  high 
accuracy  [8],  Stepwise  regression  analysis  is  a  well-known  process 
in  statistical  analysis  of  experimental  data  [7],  A  short  overview 
on  the  method  is  given  in  the  next  paragraphs.  For  more  detailed 
background  information  we  refer  to  Refs.  [7-11], 

First  of  all,  a  linear  model  is  taken  with  terms  of  all  input  vari¬ 
ables,  transformation  of  these  ones  (mostly  squared,  logarithmic, 
exponential,  and  reciprocal),  and  interactions  between  the  input 
variables.  After  the  parameter  estimation,  a  null  hypothesis  is  tested 
for  each  estimated  parameter  |,  as  described  in  Refs.  [10,12,13],  The 
null  hypothesis  is  denied  if  the  probability  of  the  test  0  is  smaller 
than  the  chosen  significance  level  0out  and  therefore,  the  tested 
coefficient  §  remains  in  the  model.  If  the  probability  is  greater  than 
the  significance  level  (0  >  0Out)  the  tested  coefficient  is  assumed 
to  be  zero  and  the  corresponding  term  is  taken  out  of  the  model. 
To  improve  numerical  stability  and  to  speed  up  regression  conver¬ 
gence,  a  second,  smaller  significance  level  0in  can  be  defined.  If  the 


Since  the  parameters  for  the  minimum  J  have  to  be  found,  Eq.  (6) 
has  to  be  minimized  by  differentiation  and  must  be  zero.  Finally,  the 
optimal  parameters  can  be  calculated  by  the  following  equation: 

§  =  (X'X)-1X'Y  (7) 

The  next  step  is  to  determine  the  standard  deviation  sd  for  each 
parameter  (with  parameter  index  6  =  1 . .  .0),  which  is  defined  by  the 
product  of  the  standard  deviation  of  the  regression  sregresSi0n  and  the 
variance  of  each  variable,  which  can  be  found  in  the  main  diagonal 
elements  of  the  covariance  matrix  as  described  in  the  following 
equation: 


s0  —  degression 


^((x'xr1)^ 


(8) 


For  the  calculation  of  sregression  (Eq.  (9)),  is  the  degree  of  freedom 
defined  by  the  difference  between  the  number  of  experimental 
0  data  and  the  number  of  parameters  0,  is  important: 


EliM-Ydf 


(9) 
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The  ratio  of  parameter  and  standard  deviation  is  the  Student- 
factor  t  of  the  null  hypothesis  test: 


By  knowing  this  value,  the  probability  can  be  calculated  by  solv¬ 
ing  the  integral  of  the  t-distribution  function: 


%  = 


rw  +  W) 

sf&nrwp) 


with  r  the  gamma  function  (Eq.  (12)): 


™-jT 


e  V  1  ds 


(11) 


(12) 


Thus,  it  is  possible  to  calculate  Vg  for  every  variable  of  the 
polynomial.  To  evaluate  the  quality  of  the  resulting  model,  a  good 
criterion  is  the  degrees  of  freedom  adjusted  R-square  (R^ ).  There¬ 
fore,  the  sum  of  squares  due  to  error  ( SSE )  has  to  be  calculated  as 
given  in  the  following  equation: 


SSE  =  ^(Yd-Yd) 


(13) 


With  the  mean  value  Yd  the  sum  of  squares  about  the  mean  (SSM) 
can  be  calculated: 


SSM=^(Yd-Yd) 


(14) 


Then,  a  measure  for  model  quality  is  developed,  the  so-called  R- 
square  (R2): 


R  =  1  "  SSM  (15) 

The  disadvantage  of  this  measure  is  that  it  always  becomes 
larger  while  adding  more  parameters  to  the  model.  This  problem 
can  be  solved  by  introducing  the  degrees  of  freedom  adjusted  R- 
square(R2dj)[10]: 


*adj  =  l- 


SSE(<£-1) 

SSM(#) 


(16) 


The  R2dj.  declines  if  additional  model  parameters  do  not  significantly 
decrease  SSE. 


3.  Experimental  work 

3.1.  Direct  steam  injection  (DSI) 

Our  team  showed  recently  [3]  that  the  investigation  of  dynamic 
behaviour  by  step  response  and  identification  of  the  corresponding 
transfer  function  is  a  promising  way. 

For  dynamic  analysis,  the  investigated  system  has  to  be  able  to 
manage  step  responses  fast  and  without  a  delay.  This  assumption 
is  not  achieved  by  state  of  the  art  in  humidification,  the  so-called 
bubbler-humidifiers.  A  typical  bubbler-humidifier  and  the  corre¬ 
sponding  working-principle  are  depicted  in  Fig.  1 . 

On  its  way  to  the  stack,  the  gas  has  to  pass  a  packed  bed  and 
absorbs  thereby  down  streaming  water.  The  variable  of  interest,  the 
relative  humidity  rh,  can  be  controlled  by  adjusting  the  ratio  of  tem¬ 
perature  in  the  humidifier  and  temperature  in  the  stack.  Since  the 
gas  leaves  the  humidifier  with  the  dew  point  temperature,  the  gas 
will  always  have  a  relative  humidity  of  rh  =  1.00  when  it  is  entering 
the  stack.  Because  of  its  low  heat  capacity,  the  gas  will  be  match¬ 
ing  to  the  stack  temperature  very  fast  in  a  heat  exchanger  and  the 
relative  humidity  will  change.  If  both  temperatures  have  the  same 
value,  the  relative  humidity  of  1.00  will  be  persist,  however  if  the 
temperature  of  the  gas  is  cooled  down  in  the  heat  exchanger  for 
example,  from  75  to  59  °C,  the  relative  humidity  will  just  decrease 
to  rh  =  0.50. 

The  disadvantage  of  this  robust  and  long-stable  system  is,  as 
mentioned  above,  the  missing  ability  to  respond  to  variations 
in  humidification  fast  enough  because  of  its  high  inert  thermal 
behaviour  and  its  large  volume.  In  order  to  identify  transfer  func¬ 
tions  from  experimental  data,  a  system  must  be  able  to  apply  ideal 
steps  with  maximum  time  delays  in  the  range  of  a  few  seconds.  A 
bubbler-humidifier  takes  for  the  mentioned  step  in  humidity  about 
4  min  and  is  therefore  inapplicable  for  dynamic  investigations. 

The  problem  is  solved  by  the  development  of  a  more  sophisti¬ 
cated  system,  enabling  the  direct  injection  of  steam  (DSI)  instead 
of  the  time-intensive  accumulation  of  water  in  a  gas.  A  draft 
of  the  experimental  set-up  used  for  all  tests  can  be  found  in 
Fig.  2. 

The  principle  of  DSI  is  the  connection  of  a  steam  generator  at  the 
gas  supply  and  the  controlling  of  the  entering  steam  amount  by  a 
mass  flow  controller  (F1C 1  in  Fig.  2 ),  which  will  be  explained  later,  in 
front  of  the  mixing  point  between  steam  and  dry  gas.  Heating  tapes 
keep  the  temperature  constant  and  so  no  water  can  condense  till  the 
temperature  of  the  wet  gas  is  matched  with  the  stack  temperature 
in  a  heat  exchanger.  The  required  amount  of  steam  is  calculated  by 
the  test  bench  software  in  real-time  using  Eqs.  (1)— (4). 


Fig.  2. 


iiagram  of  humidification  by  DSI. 
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As  shown,  it  is  possible  to  calculate  the  required  steam  amount 
only  by  measuring  the  mass  flow  of  dry  gas  in  front  of  the  mix¬ 
ing  point  (FIC  2  in  Fig.  2),  the  pressure  of  the  gas  at  stack  inlet, 
and  the  temperature  (TIC  2  in  Fig.  2)  at  stack  inlet.  For  tempera¬ 
ture  measurement,  coolant  temperature  is  more  suitable  than  gas 
temperature.  The  thermal  energy  of  the  coolant  is  2  orders  of  mag¬ 
nitude  higher  than  the  thermal  energy  of  the  reactants.  Due  to  a 
large  available  heat  transfer  area  and  the  good  thermal  conductiv¬ 
ity  of  the  used  materials,  the  gas  temperature  is  rapidly  adjusted  to 
coolant  temperature  at  the  stack  inlet. 

Thermocouples  and  pressure  sensors  are  cost-effective  and 
accurate  measuring  sensors  and  were  therefore  chosen.  Measure¬ 
ments  of  the  mass  flow  are  more  complex  and  were  realized  by 
coriolis  mass  flow  controller  (CMFC)  CORI-Flow  from  Bronkhorst. 
These  sensors  measure  the  created  phase  shifting  of  two  oscillat¬ 
ing  curved  tubes  by  coriolis  force.  When  a  mass  flow  is  passing  the 
tubes,  CMFCs  have  great  advantages  like:  real-time  measurement, 
independence  of  material  properties  like  density  or  viscosity  and 
a  high  accuracy  with  a  error  of  just  about  0.2%.  A  more  detailed 
description  of  the  coriolis  measuring  principle  is  given  in  Ref.  [14], 
For  the  validation  of  the  proposed  humidification  by  direct  steam 
injection  precise  measurement  of  gas  humidity  is  necessary.  There¬ 
fore,  two  common  used  capacitive  hygrometers  HTS  32D  from 
Rotronic  with  a  specified  accuracy  of  2%  in  relative  humidity  and  a 
chilled  mirror  dew  point  hygrometer  (CMDPH)  DP3-D  from  MBW 
with  an  accuracy  of  less  than  0.1%  relative  humidity  are  used. 

CMDPHs  are  using  the  principle  of  cooling  down  a  mirror  ther- 
moelectrically  by  a  peltier  device  until  the  mirror  reaches  the  dew 
point  of  the  surrounding  gas.  At  this  temperature,  the  water  will 
begin  to  form  condensate  on  the  mirror  resulting  in  an  improved 
absorption  of  emitting  light.  The  decreased  light  reflection  is  notic¬ 
ing  by  a  detector  and  a  software  is  able  to  calculate  the  relative 
humidify  using  Eqs.  (1)  and  (2).  The  used  CMDPH  is  accurately 
calibrated  by  the  manufacturer  and  is  used  as  a  reference  during 
validation  of  the  DSI  method. 

3.2.  Cathode  humidity  cycling 

All  fuel  cell  tests  were  done  on  a  test  bench  of  Daimler 
using  a  PEMFC  stack  with  metallic  bipolar  plates  and  commer¬ 
cial  Membrane-Electrode-Assembly  (MEAs).  The  stack  design  was 
developed  by  Daimler  and  was  already  used  in  previous  work  [3], 
because  of  its  strong  similarity  to  stacks  used  in  Daimler’s  prototype 
fuel  cell  car  F  600. 

The  stack  consisted  of  ten  cells,  with  an  active  area  of  300  cm2 
each  and  is  integrated  under  a  pressure  of  3  bar.  All  experiments 
are  carried  out  in  a  fully  software  controlled  Daimler’s  Fuel-Cell- 
Stack-Test-Bench.  The  cells  are  operated  with  cleaned  air  as  cathode 
gas  and  a  hydrogen/nitrogen  mixture  as  anode  gas.  All  gas  flows 
are  controlled  by  two  mass  flow  controllers,  one  for  small  mass 
flows  and  one  for  large  mass  flows,  to  ensure  a  high  accuracy  in  the 
whole  operating  area.  The  pressures  of  both  gases  are  adjusted  by 
digital  controlled  valves  at  the  outlet  and  the  stack  temperature  is 
controlled  by  varying  the  mass  flow  of  cooling  water. 

The  test  bench  allows  measurements  in  the  major  operating  area 
of  fuel  cell  powered  drive  trains,  in  detail  this  are:  a  gauge  pres¬ 
sure  p  until  2.6  bar,  a  stoichiometry  X  (ratio  of  supplied  reactant  to 
chemically  needed  reactant)  from  1.4  to  3,  a  current  density  i  in  the 
range  of  0.4-1.6Acm-2,  a  stack  inlet  temperature  Tin  between  40 
and  80  °C,  and  a  difference  between  inlet  and  outlet  temperature 
AT  of  0-10  K.  Restrictions  had  just  been  made  in  regard  to  comply 
with  a  maximum  pressure  difference  of  500  mbar  between  anode 
and  cathode.  A  further  input  parameter  in  all  experiments  is  the 
fraction  of  H2  in  the  anode  gas  Ch2,  since  the  diffusion  of  air  from 
cathode  to  anode  is  a  typical  phenomenon  of  PEMFCs.  In  automotive 


Table  1 

Experimentally  investigated  conditions 
Input  <9 

Current  density  i 

Cathode  gas  stoichiometry  Xc 

Anode  gas  stoichiometry 
Cathode  gas  pressure  pc 

Anode  gas  pressure  pa 

Relative  humidity  anode  gas  rha 

Inlet  temperature  Tin 

Outlet-inlet  temperature  AT 

Hydrogen  concentration  ch2 


0.33-1 

40-75 

0-18 

0.4-1 


drive  trains  this  accumulation  of  N2  is  intensified  by  the  recycling 
of  anode  gas  by  a  H2-loop  and  the  fraction  of  hydrogen  in  the  anode 
Ch2  is  therefore  a  point  of  interest.  In  Table  1,  the  experimentally 
investigated  operation  conditions  are  summarized. 

For  covering  the  listed  drive  train  relevant  operating  area 
reasonably,  101  operating  conditions  are  investigated  without 
applying  harmful  conditions  to  the  stack. 

A  three  level  humidity  cycle  (rhc  =  0.33,  0.67,  and  1.00)  is  con¬ 
ducted  as  described  in  Fig.  3  to  stimulate  cathode  gas  relative 
humidity-related  voltage  dynamics.  With  a  sample  rate  of  1  Hz,  the 
data  are  recorded,  exported  to  MATLAB,  and  analysed  in  the  next 
chapter. 

4.  Experimental  results 

4.1.  DSI  verification 

The  sensors  for  the  measurements  of  relative  humidity  work 
like  expected  very  slowly,  about  1  h  must  be  scheduled  for  getting 
an  exact  value.  Therefore,  the  direct  measurement  of  gas  humidity 
is  not  possible.  In  Fig.  4,  the  measured  relative  humidity  is  plot¬ 
ted  against  the  mass  flow  of  steam  injected  into  dry  air  for  an  air 
flow  of  100  N1  min-1 .  Furthermore,  the  calculated  relative  humidity 
(Eqs.  (1 )— (4))  using  the  measured  steam  mass  flow,  air  mass  flow, 
temperature,  and  pressure  is  given  in  Fig.  4. 

The  capacitive  sensors  are  only  used  for  having  a  comparison 
and,  are  because  of  a  high  average  deviation  of  about  5%  relative 
humidity  and  a  maximum  deviation  of  30%  points,  are  not  suffi¬ 
cient  for  precise  measurements  and  thus  not  used  for  calibration, 
whereas  the  CMDPH  had  only  an  average  deviation  of  0.5%  and 
the  deviations  are  low  in  the  whole  measuring  range.  For  demon¬ 
strating  the  high  accuracy  of  the  DSI  method  in  a  large  operating 
area,  the  CMDPH  was  tested  also  for  smaller  mass  flows  and  the 
results  are  plotted  in  Fig.  5.  The  lines  in  the  upper  part  of  this  fig¬ 
ure  show  the  experimental  and  calculated  data  for  a  mass  flow  of 


Fig.  3.  Realized  steps  in  humidity  for  every  run. 
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Fig.  4.  Comparison  of  measured  and  calculated  relative  humidity  for  a  mass  flow  of 
100  N1  min-1  dry  air. 


20Nlmin_1  of  dry  air  and  the  ones  in  the  lower  part  for  a  mass 
flow  of  50  N1  min-1  of  dry  air.  The  average  deviations  are  again  just 
low  and  as  a  consequence  the  CMDPH  was  used  for  the  calibration 
before  every  measurement. 

A  very  accurate  and  most  of  all,  very  dynamic  method  for  vary¬ 
ing  the  gas  humidity  is  available  now.  Impressive  proved  by  a  time 
duration  of  just  2  s  to  manage  the  step  from  0.50  to  1.00  in  relative 
humidity  and  thus  a  strong  enhancement  to  previous  methods.  Ref. 
[9]  calls  a  system  suitable  for  dynamic  investigations  if  it  is  2-10 
times  faster  than  the  step  response.  In  our  case,  the  voltage  needs 
at  least  30  s  or  15  times  of  the  time  of  the  system,  so  this  rule  of 
thumb  is  fulfilled  and  dynamic  analysis  is  possible. 

4.2.  Cell  voltage  response 

Similar  to  Ref.  [15],  the  experimental  data  of  the  stack  perfor¬ 
mance  shows  the  effect  that  the  water  content  of  the  membrane 
is  more  related  to  the  cathode  humidification  than  to  the  anode 
humidification.  That  means  nearly  no  stationary  voltage  effect  of 
anode  humidification  is  observable  even  at  low  cathode  gas  humid¬ 
ification.  Besides,  no  dynamic  voltage  can  be  seen  for  steps  in  the 
anode  humidity  and  for  this  reason,  just  the  cathode  is  investigated 
in  the  further  work. 

In  different  operation  points  (exemplary  in  figure  caption 
Figs.  6  and  7)  of  the  PEMFC  system  measurements  were  accom¬ 
plished  with  the  same  cathode  humidity  cycle  (Fig.  3).  Afterwards 
the  dynamic  in  the  mean  voltage  is  analysed  according  to  cathode 
humidity. 

Therefore,  a  correct  predicting  model  had  to  be  able  to  consider 
not  only  the  static  behaviour,  but  also  the  dynamic  behaviour  of  the 
stack  voltage  in  dependence  on  operating  parameters. 


Fig.  5.  Measured  and  calculated  relative  humidify  for  20  and  50Nlmin_1  flow  of 


Fig.  6.  Variation  of  humidity  and  measured  step  response  of  voltage  in  case  of  an 
operating  point  with  strong  dynamic  (i  =  1.4Acm-2,  2.c  =  1.9,  2.a  =  1.8,  pc  =  1.8bar, 
Pa  =  2.2  bar,  rha  =  0.33,  T=  55  °C,  AT=  10  K,  ch2  =  1 ). 


One  of  the  most  interesting  results  of  the  experimental  data 
analysis  is  the  fact  that  two  completely  different  behaviours  of  the 
system  are  observed  when  the  cathode  humidification  is  varied.  In 
Fig.  6,  the  mean  voltage  of  the  ten  cells  shows  a  strong  correlation 
to  the  humidity  of  the  cathode,  however,  in  Fig.  7,  nearly  no  influ¬ 
ence  of  the  relative  humidity  is  observable.  Also  important  is  the 
fact  that  the  voltage  in  Fig.  6  shows  a  significant  dynamic  behaviour 
after  each  step.  Fig.  6  shows  that  the  cathode  humidity  of  this  oper¬ 
ation  point  has  a  significant  influence  in  the  static  and  dynamic 
behaviour  of  the  mean  voltage.  Flowever,  Fig.  7  shows  an  operation 
point  with  almost  no  influence  of  humidity  in  cell  voltage.  In  depen¬ 
dence  of  the  step  height,  it  took  5-33  min  till  the  new  steady  state 
is  reached.  In  other  publications  the  same  behaviour  of  the  cath¬ 
ode  humidity  according  to  the  different  amplitude  of  the  mean  cell 
voltage  in  different  operation  points  can  be  observed  [16-18], 

5.  Modelling  of  dynamic  stack  voltage 

The  total  of  606  steps  are  divided  into  a  set  of  174  steps  for 
parameter  identification  and  the  remaining  ones  for  the  indepen¬ 
dent  validation  of  the  gained  model.  For  every  of  these  174  steps,  the 
time  constants  were  identified  by  an  optimization  function  using 
the  least  square  method. 

The  focus  in  this  work  is  set  on  the  modelling  of  the  dynamic 
portion  of  cell  voltage.  This  procedure  saves  the  already  done  mod¬ 
elling  of  the  static  behaviour  and  allows  the  fusion  of  both  portions 
at  the  end.  The  dynamic  voltage  l/dynamic.  defined  as  the  differ¬ 
ence  between  measured  voltage  IWasured  and  steady  state  voltage 


Fig.  7.  Variation  of  humidity  and  measured  step  response  of  voltage  in  case  of 
an  operating  point  with  low  dynamic  (i  =  0.9Acm“2,  2,c=3,  1-a  =  3,  pc  =  1.5bar, 
pa  =  1.5  bar,  rha  =  0.33,  T= 70  °C,  F=  10  K,  cH2  =  1). 
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Fig.  8.  Comparison  of  measured  dynamic  voltage  and  dynamic  voltage,  simulated 
by  a  transfer  function  with  TD  =  633.9  and  Ka — 22.2  (i  - 1.4  A  cm-2,  Xc  =  3.5,  X3  -  1.5, 
pc  =  2.1  bar,  pa  =  1.7  bar,  rha = 33%,  T=  70  °C,  T=  0  K,  Ch2  =  1 ,  rhc>1  =  1.00,  rhc,2  =  0.33). 


Fig.  9.  identified  vs.  estimated  KD  using  Eq.  (7)  (x  =  fitted  KD,  dotted  line  =  ideal 
regression). 


^static  (Eq.  (17)),  was  identified  for  every  run  and  analysed  in  view 
of  describing  its  behaviour: 


^dynamic  —  ^measured 


-  ^static 


(17) 


The  dynamic  voltage  shows  in  all  cases  the  typical  behaviour  of 
a  lst-order  proportional-differential  transfer  function  PDT,  which 
is  described  by  the  following  equation: 


C(s)  = 


Kps 

fm+i 


(18) 


The  dynamic  voltage  was  simulated  by  using  this  transfer  func¬ 
tion.  In  this  respect,  /<□  is  characterising  the  height  of  the  step 
response  and  TD  the  time  delay  till  the  steady  state  is  reached 
again.  One  of  these  simulated  step  responses  in  comparison  to  the 
measured  dynamic  behaviours  is  shown  in  Fig.  8. 

The  simulated  step  in  cathode  gas  humidity  by  an  identified 
PDT]  transfer  function  fits  very  well  the  measured  dynamic  voltage 
and  was  therefore  a  promising  way  to  characterize  the  fuel  cell. 


5.1.  Linear  modelling 


First,  a  simple  linear  dynamic  model  is  investigated.  To  describe 
the  dynamic  cell  voltage,  the  previously  introduced  PDT]  transfer 
function  (Eq.  (18))  with  constant  coefficients  Ko  and  Tp  is  used. 
The  174  steps  for  parameter  estimation  can  be  fitted  best  with  a 
Kp  =  —9.1  and  a  TD  =  633.9. 


5.2.  Nonlinear  modelling 

The  experimental  analysis  of  all  174  steps  individually,  a  strong 
correlation  between  KD  and  Tp  can  be  observed.  Due  to  this  high 
correlation,  the  delay  time  is  fixed  to  the  value  of  the  linear  model 
as  Tp  =  633.9.  The  estimated  values  for  KD  vary  in  a  range  from 
0  to  -60.  However,  the  step  height  Kp  depends  nonlinear  on  the 
input  variables.  To  develop  a  precise  and  simple  polynomial,  witch 
describes  the  dependency  of  KD  regarding  the  operating  conditions, 
a  stepwise  multiple-regression  analysis  is  used. 

In  the  first  step,  the  matrix  X  is  set  up.  In  the  first  column,  ones 
are  placed  to  include  a  constant  term  in  the  polynomial.  Columns 
2-10  are  filled  with  the  untransformed  inputs,  as  they  are  i,  Xc, 
X3,  pc.  Pa,  rha,  Tin,  AT,  and  Ch2.  Columns  11-19,  with  the  logarithm 
of  the  inputs,  columns  20-28,  with  the  exponential  transforma¬ 
tion,  in  columns  29-37,  the  reciprocal  inputs,  and  columns  38-46, 
the  square  root  of  the  9  inputs.  Secondly,  the  matrix  X  is  expanded 
in  the  columns  47-1036  with  all  two-way-interactions  of  columns 
2-46  like  iXc,  iX3,  ipc,  etc.  Finally,  columns  1037-1081  are  added 
with  the  quadratic  values  of  columns  2-46  to  gain  a  full  quadratic 
polynomial. 


In  the  second  step,  the  stepwise  multiple-regression  analysis 
is  initialized  considering  |i-46  to  be  nonzero  and  therefore  be 
included  terms  in  the  polynomial.  All  other  coefficients  £47-1081  are 
initialized  to  be  zero  and  therefore  off  the  polynomial.  The  remove 
threshold  '/'out  is  set  to  '/'out  =  0.05  and  the  enter  threshold  •Pin  is 
set  to  ipout  =  0.1. 

Thirdly,  stepwise  multiple-regression  is  started.  The  probabil¬ 
ity  of  the  null  hypothesis  '/'g  of  all  coefficients  is  calculated 
as  described  in  chapter  2  using  Eqs.  (5)-(12).  After  that,  all  non¬ 
significant  terms  are  removed  and  significant  terms  which  have 
been  previously  off  are  included  to  the  polynomial. 

Finally,  the  third  step  is  repeated  12  times  while  no  change  in 
individual  term  significance  can  be  observed. 

The  result  of  the  stepwise  multiple-regression  analysis,  using 
the  untransformed  9  inputs  and  the  proposed  transformation,  are 
12  significant  terms  as  they  are  summarized  in  Table  2. 

The  Rjd j  of  the  stepwise  multiple-regression  is  quite  good  with 
f?adj  =  0.9849  which  can  be  visualized  with  a  estimated  XD  (calcu¬ 
lated  with  the  model  shown  in  Table  2)  vs.  identified  Xq  (identified 
with  the  measured  data)  plot  as  drawn  in  Fig.  9  where  every  cross 
represents  a  different  data  set  in  an  own  operation  point  of  the 
PEMFC  system. 


6.  Model  validation 

Always  if  experimental  based  system  identification  methods  are 
used,  the  developed  models  must  be  validated,  using  independent 
experiments.  From  the  606  experimental  step  data  sets,  just  174 
have  been  used  during  parameter  estimation  process.  To  do  the 
validation  of  the  postulated  polynomial  to  calculate  KD  as  a  function 
of  PEMFC  operation  conditions,  the  432  not  jet  used  experimental 
data  sets  are  used. 


Table  2 

Result  of  stepwise  multiple-regression  analysis  of  PEMFC  operating  conditions  influ¬ 
ence  at  Ko 
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time  (s) 


Fig.  10.  Comparison  of  measured  and  simulated  dynamic  voltage  of  a  run  with 
low  dynamic  voltage  (i  =  1 A  cm-2,  Xc  =  1.9,  T.a  =  1.4,  pc  =  1.4  bar,pa  =  1.6  bar,  rha  =  0.50, 
r=55°C,r=10K,CH2  =  1). 


time  (s) 


Fig.  11.  Comparison  of  measured  and  simulated  dynamic  voltage  of  a  run  with  high 
dynamic  voltage  (i  =  0.9  A  cm-2,  Xc  =  1.4,  T.a  =  1.4,  pc  =  1.8  bar,  pa  =  1.8  bar,  rha  =  0.66, 
T=55°C,  AT=10K,Ch2  =1). 

A  calculated  degrees  of  freedom  adjusted  R-square  for  the  vali¬ 
dation  data  of  R^d|  =  0.9689  indicates  a  good  generalization  of  the 
developed  polynomial.  The  match  between  measured  and  calcu¬ 
lated  dynamic  voltage  using  Kb  from  the  polynomial  in  Eq.  (18) 
for  two  representative  runs  is  given  in  Figs.  10  and  11.  In  contrast 
to  Figs.  6  and  7  shows  a  pure  dynamic  of  the  mean  cell  voltage  in 
different  operation  points.  Fig.  10  has  a  lower  dynamic  with  about 
45  mVpp  (mV  peak  to  peak)  in  comparison  to  Fig.  11  with  a  higher 
dynamic  with  about  80  mVpp. 

It  is  obviously,  that  the  introduced  method,  the  combination 
of  transfer  function  and  nonlinear  regression  of  Kp,  describes  the 
dynamic  voltage  very  satisfying.  Especially,  the  fact  that  different 
dynamic  behaviours,  the  second  run  curve  has  a  dynamic  volt¬ 
age  twice  as  big  as  the  first  one,  were  simulated  correctly  only  by 
knowing  the  operating  conditions  is  a  success. 

Particularly  the  comparison  of  the  pure  static,  the  linear 
dynamic,  and  the  nonlinear  dynamic  model  give  a  good  impression 
of  the  realized  improvement  in  PEMFC  modelling.  For  comparing 
these  three  models  and  thus  evaluate  our  new  method,  the  error 
between  calculated  dynamical  voltage  and  simulated  one  is  anal¬ 
ysed  statistically.  The  residual  histograms  in  Figs.  12-14  are  similar 
to  the  Gaussian  Distribution  and  the  shape  of  the  Distribution  deliv- 


residual 


Fig.  12.  Residual  histogram  of  the  static  model. 


Fig.  13.  Residual  histogram  of  the  linear  model. 


residualnonljii 


Fig.  14.  Residual  histogram  of  the  nonlinear  model. 


ers  an  estimation  of  the  simulation  quality.  The  first  plot  (Fig.  12),  is 
performed  without  a  consideration  of  the  dynamic  behaviour  and 
is  therefore  a  good  reference  for  the  dynamic  simulations  of  stack 
voltage. 

Figs.  13  and  14  show  the  mean  error  ME  could  be  minimized  at 
about  35%  by  using  a  PDTi  -element  with  a  constant  KD  and  a  further 
significant  improvement  of  about  69%  is  achieved  by  replacing  this 


Table  3 

Results  of  model  validation 


Model  ME  (mV)  1  “  (%)  (mV)  i  _  (%)  1  -  ■J/’41  (%) 

Static  1.67  0.0  8.28  0.0  0.8943  0 

Linear  1.09  34.7  5.55  33.0  0.9353  38.8 

Nonlinear  0.34  79.6  4.27  48.4  0.9689  70.6 
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time  (s) 


Fig.  15.  Comparison  of  simulated  and  measured  mean  stack  voltage  of  the  nonlin¬ 
ear  model  (i  =  1  A  cm-2,  kc  =  1.9,  =  1.4,  pc  =  1.4 bar,  pa  =  1.6bar,  rha  =  0.50,  T=55  °C, 

AT  =  10  K,  Ch2  =1). 

constant  I(D  by  using  the  polynomial  given  in  Table  2.  Compared, 
the  pure  static  model  a  total  reduction  of  mean  error  of  80%  can  be 
realized.  A  final  mean  error  of  nearly  zero  (ME  =  0.34  mV)  proves, 
how  accurate  the  dynamic  voltage  can  be  simulated.  The  reduction 
of  standard  error  from  8.28  to  4.27  mV,  if  the  nonlinear  model  is 
added  to  the  static  model,  equals  nearly  a  bisection  of  the  error. 

In  Table  3,  the  results  of  the  model  validation  are  summarized. 
Nearly  90%  of  overall  precision  is  related  to  the  stationary  model  as 
the  statjc  value  indicates.  38.8%  of  gas  humidity-related  voltage 
dynamics  can  be  described  with  a  simple  PDTi  transfer  function 
using  constant  coefficients.  A  highly  precise  nonlinear  model  is 
developed,  which  reduces  model  errors  by  70.6%.  Such  a  model  is 
able  to  fulfill  highest  demands  for  system  simulation  and  power 
train  optimization  purposes. 

In  Fig.  15,  the  mean  cell  voltage  of  the  investigated  PEMFC  stack 
is  plotted  for  a  complete  cathode  humidity  cycle.  The  good  perfor¬ 
mance,  indicated  by  the  residual  histogram,  can  be  visualized.  The 
noise  of  the  simulated  voltage  is  related  the  noise  of  the  measured 
9  inputs,  which  are  fed  directly  into  the  model.  If  stable  inputs  are 
used,  the  simulated  voltage  is  smooth. 

7.  Summary  and  outlook 

A  general  method  for  simulating  the  gas  humidity-related  volt¬ 
age  dynamics  of  a  PEMFC  stack  is  developed  and  introduced.  For 
this  purpose,  a  new  humidification  system  has  been  developed  and 
validated  against  highly  precise  CMDPHs.  The  DSI  allows  dynamic 
investigations  of  a  PEMFC  without  the  need  of  humidity  mea¬ 
surements  and  is  able  to  provide  experimental  data  to  identify 
a  dynamic  model.  Step  responses  of  the  mean  cell  voltage,  after 


varying  the  humidification  of  the  cathode,  are  measured  and  the 
resulting  behaviour  is  well  described  by  a  PDTI  -element.  The  pro¬ 
portional  constants  Kd.  gained  from  parameter  estimation,  are 
used  in  combination  with  the  operating  conditions  to  identify  a 
2nd-order  polynomial  by  a  stepwise  multiple-regression  analy¬ 
sis,  which  approximates  accurately  the  dynamical  behaviour  of  all 
steps.  Transfer  function  and  regression  equation  are  successfully 
validated.  The  fit  between  simulated  and  experimental  data  is  high 
and  thus  a  novel  dynamic  model  for  the  stack  voltage  is  successfully 
introduced. 

Possible  applications  of  this  new  method  are,  for  example,  the 
substitution  of  experiments  by  simulations,  and  the  development 
of  an  operating  strategy  of  PEMFC  powered  drive  trains. 

Further  work  should  be  done,  by  investigate  more  operating 
conditions  to  improve  the  valid  range  of  the  model. 
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